Profr. Julio César Ramírez Pacheco

El concepto de entropía

La entropía se puede interpretar como un indicativo de la complejidad en una serie de datos aleatorios \(\{X_1, X_2, \ldots, X_N\}\). La entropía usualmente se calcula utilizando una función de masa de probabilidad \(p_j\) que cumple con las siguientes propiedades:

  • \(p_j \ge 0, \, j \in 0,1, \ldots, N\)
  • \(\sum_{j=1}^N{p_j} = 1\)

R tiene múltiples paquetes y funcionalidades que permiten estimar la pmf de un conjunto de datos como el descrito anteriormente. El histograma es una herramienta que permite estimar la pmf de un conjunto de datos. En el ejemplo que sigue se muestra la forma de estimar la pmf en R:

set.seed(1234)                       # Para hacer el análisis reproducible
datos     <- rnorm(512,0,1)          # Se generan 512 valores normales
histogram <- hist(datos, plot=FALSE) # Se calcula el histograma
pmf       <- histogram$counts/sum(histogram$counts)  # Se calcula la pmf
sum(pmf)                             # Se verifica que cumpla con las propiedades
## [1] 1

Ejercicios:

  • Estimar la pmf utilizando el utilizando los paquetes ASH y KernSmoooth.

R=

Usando el paquete ASH, la pmf se calcula de la siguiente forma.

pacman::p_load(ash)

# Se calcula el histograma desplazado promediado univariado (en inglés, univariate averaged shifted histogram o ASH) usando un kernel polinomial, que es el que usa por defecto el paquete.
histogram_ash <- ash1(bin1(datos))
## [1] "ash estimate nonzero outside interval ab"
pmf <- histogram_ash$y / sum(histogram_ash$y)
sum(pmf)
## [1] 1

Usando el paquete KernSmooth, la pmf se calcula de la siguiente forma.

pacman::p_load(KernSmooth)

# Estimación de densidad de kernel agrupada (Binned Kernel Density Estimate) sobre los datos
bkde_kernsmooth <- bkde(datos)

pmf <- bkde_kernsmooth$y / sum(bkde_kernsmooth$y)
sum(pmf)
## [1] 1
  • ¿Cuál es la ventaja de utilizar los métodos anteriores sobre el histograma?

R=

Los métodos que ofrecen estos paquetes son mucho más robustos que el histograma para calcular la pmf. Además, permiten aproximar la función de distintas formas, ya sea usando distintos kernels o llevando a cabo la estimación para una función en 2D.

  • Utilizando el comando hist y los paquetes ASH y KernSmooth verifique el tiempo requerido para estimar la densidad de una serie de datos Gaussianos con \(\mu=1\) y varianza \(\sigma^2=1\) y longitudes \(N=2^i, \, i=8,9,10, 11, \ldots 16.\) (es necesario incluir un gráfico en highcharter)

R=

Para calcular el tiempo de ejecución, se usa la función sys.time en cada ejecución de las 3 funciones con los distintos valores de \(i\).

Los tiempos de ejecución se grafican con el paquete highcharter.

pacman::p_load(highcharter)

# Graficar los resultados usango highcharter
x <- 8:16
highchart() %>%
  hc_add_series(cbind(x, res_hist), name = "Función hist") %>% 
  hc_add_series(cbind(x, res_ash1), name = "Función ash1") %>%
  hc_add_series(cbind(x, res_bkde), name = "Función bkde") %>%
  hc_add_theme(hc_theme_smpl()) %>%
  hc_title(text = "Tiempo de ejecución") %>%
  hc_subtitle(text = "IT0322 - Teoría de la información") %>%
  hc_xAxis(title = list(text = "i")) %>%
  hc_yAxis(title = list(text = "Tiempo"))

La entropía utilizando el histograma

Volviendo de nuevo al ejemplo anterior, podemos estimar la entropía de Shannon, utilizando la pmf obtenida mediante el histograma y así obtener un estimador empírico de la entropía de Shannon. A continuación mostramos la forma de obtener la entropía de una serie de datos obtenida en ventanas independientes o contiguas de longitud \(512\):

set.seed(1234)
datos        <- rnorm(32768)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Preguntas

  1. ¿Porqué existen valores discontinuos en la entropía?

R=

Porque existen puntos cuyos valores de pmf son 0, por lo que al calcular el logaritmo de la entropía sobre estos valores provoca que el valor resultante esté indefinido (en R esto sería obtener un NaN), provocando que estos valores se vean como discontinuidades en la gráfica.

  1. ¿Con qué código soluciona el problema de las discontinuidades?

R=

Si al momento de calcular la suma de la entropía usando la función sum excluyéramos estos valores con el parámetro na.rm = TRUE.

  1. Ahora calcule la entropía de una serie de datos normales, denotados por \(X_t\) (con \(\mu=0\) y \(\sigma=2\)) pero ahora añádanle (súmenle) una segunda función \(r_t\), es decir, hallen la entropía de la serie \(Y_t = X_t+r_t\) con \(r_t\) definida por: \[ r_t = \begin{cases} \sigma/4 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 4, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

  1. Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma/2 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 2, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

  1. Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

  1. ¿Tiene algún efecto la longitud del salto en \(r_t\) en la forma de la entropía? Explique.

R=

  1. ¿Qué sucede ahora si \(r_t\) es de la forma: \[ r_t = \begin{cases} \sigma & 15872 \le t\le 16896\\ 0 & \mbox{otro caso} \end{cases} \] ?
# Funcion r_t
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

  1. Repita los pasos 3-7 pero ahora usando la entropía de Harvda con parámetro \(\alpha=3\) y \(\alpha=9\). ¿Qué efecto tiene \(\alpha\)?

Repitiendo el paso 3, con ambos valores de alpha para la entropía de Havrda se obtiene lo siguiente:

Usando \(\alpha = 3\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 4, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 3)) * log(sum(pmf^3) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Usando \(\alpha = 9\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 4, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 9)) * log(sum(pmf^9) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Repitiendo el paso 4, con ambos valores de alpha para la entropía de Havrda se obtiene lo siguiente:

Usando \(\alpha = 3\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 2, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 3)) * log(sum(pmf^3) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Usando \(\alpha = 9\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma / 2, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 9)) * log(sum(pmf^9) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Repitiendo el paso 5, con ambos valores de alpha para la entropía de Havrda se obtiene lo siguiente:

Usando \(\alpha = 3\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 3)) * log(sum(pmf^3) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Usando \(\alpha = 9\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 9)) * log(sum(pmf^9) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

PASO 6

Repitiendo el paso 7, con ambos valores de alpha para la entropía de Havrda se obtiene lo siguiente:

Usando \(\alpha = 3\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 3)) * log(sum(pmf^3) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Usando \(\alpha = 9\):

# Funcion r_t
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)

set.seed(1234)
datos        <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  # Entropia de Havrda con alpha = 3
  entropies[i] <- -1*((1 / (1 - 9)) * log(sum(pmf^9) / sum(pmf)))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Entropía utilizando ventanas deslizantes

El cálculo de la entropía por ventanas independientes dada arriba resulta útil en casos en dónde la función no tiene dependencia en los valores futuros (descorrelacionadas). Para el caso de funciones correlacionadas, el cálculo de la entropía por ventanas deslizantes traslapadas resulta útil para descubrir ciertas fenomenologías en los datos. El cálculo por ventanas deslizantes de tamaño \(W\) se realiza sobre una secuencia de datos \(X_1, X_2, \ldots, X_N\). La ventana (\(W\le N\)) se va deslizando sobre los datos con factor \(\Delta\) y de esta forma subconjuntos de los datos \(X_i\) toman la siguiente forma: \[ X(m; W, \Delta) = x_j \times \Pi(\frac{t-m\Delta}{W}-\frac{1}{2}), \] donde \(m\Delta \le j \le m\Delta + W\) y \(m=0,1,2, \ldots\). Finalmente se puede graficar \(nW + \Delta, n=1,2,3, \ldots\) contra las entropías y verificar algún patrón en los datos.

Ejercicios

  • Implementar en R la metodología del cálculo de la entropía de Harvda normalizada por ventanas deslizantes. La función debe tener la forma harvda_deslizante(datos, w.length=512, s.factor=10, a.parameter=0.8, ent.type=c("hist", "ash", "kern")), donde w.length es la longitud de la ventana, s.factor es el factor de deslizamiento y a.parameter es el parámetro \(\alpha\) de la entropía de Harvda. Además la función puede calcular la entropia usando el histograma, por el método ash o por alguna metodología kernel (con el parámetro ent.type).
  • Aplicar la entropía calculada con los datos generados anteriormente, es decir:
    • Los pasos \(3,4,5\) y \(7\) en donde los datos se dan como \(X_t+r_t\).
    • ¿Tiene algún efecto el parámetro \(a\) en la forma de la entropía? ¿Tiene alguna ventaja calcular la entropía de Harvda por otro método diferente al histograma?
    • Además del histograma y los estimadores tipo kernel, existen otros métodos para estimar la distribución de una serie de datos. Investigue: ¿en qué consiste la entropía de permutación?